use "KK_efficacy_replication_data.dta", clear
gen gap = choose_jnj-choose_pfizer
gen group = .
replace group = 1 if vax_control == 1
replace group = 2 if efficacy_symptomatic == 1
replace group = 3 if efficacy_severe == 1
replace group = 4 if efficacy_hospitalization == 1
replace group = 5 if efficacy_combined == 1


* Figure 1 *
* Preferences *
reg choose_jnj vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined , nocons
estimates store A
reg choose_pfizer vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined , nocons
estimates store B
coefplot (A, symbol(circle)) (B, symbol(square)), xlabel(0 "0%" .1 "10%" .2 "20%" .3 "30%" .4 "40%" .5 "50%" .6 "60%" .7 "70%" .8 "80%") ///
ylabel(1 "No efficacy information" 2 "Symptomatic illness: 66% vs. 91%" 3 "Severe illness: 85% vs. 95%" 4 "Hospitalization/death: both > 95%" 5 "Efficacy combined") ///
graphregion(color(white)) ti("") legend(order(2 "Vax A" 4 "Vax B")) saving(choose_combined_all.gph, replace) ti("A", position(2)) scale(.75) xti(" " "% Choosing")

* Preference Gap *
reg gap vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined , nocons
coefplot, xlabel(.2 "20%" .1 "10%" 0 "0%" -.1 "-10%" -.2 "-20%" -.3 "-30%" -.4 "-40%" -.5 "-50%" -.6 "-60%" -.7 "-70%") ///
ylabel(1 "No efficacy information" 2 "Symptomatic illness: 66% vs. 91%" 3 "Severe illness: 85% vs. 95%" 4 "Hospitalization/death: both > 95%" 5 "Efficacy combined") ///
graphregion(color(white)) ti("") legend(off) saving(gap_all.gph, replace) scale(.75) xline(0) ti("B", position(2)) xti(" " "Preference Gap (% Vax A - % Vax B)")
oneway gap group, bonferroni tab
ttest gap if efficacy_symptomatic == 1|efficacy_severe == 1, by(efficacy_severe)
ttest gap if efficacy_symptomatic == 1|efficacy_severe == 1, by(efficacy_severe)
ttest gap if efficacy_symptomatic == 1|efficacy_combined == 1, by(efficacy_combined)
oneway gap group, bonferroni tab
* Combined figure *
graph combine "choose_combined_all.gph" "gap_all.gph", row(2) ysize(11) xsize(8.5) graphregion(color(white))

* Figure 2 *
* Combined -- Unvaccinated *
reg choose_jnj vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined if unvaccinated == 1 , nocons
estimates store A
reg choose_pfizer vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined if unvaccinated == 1, nocons
estimates store B
coefplot (A, symbol(circle)) (B, symbol(square)), xlabel(0 "0%" .1 "10%" .2 "20%" .3 "30%" .4 "40%" .5 "50%" .6 "60%" .7 "70%" .8 "80%") ///
ylabel(1 "No efficacy information" 2 "Symptomatic illness: 66% vs. 91%" 3 "Severe illness: 85% vs. 95%" 4 "Hospitalization/death: both > 95%" 5 "Efficacy combined") ///
graphregion(color(white)) ti("")  saving(choose_combined_unvax.gph, replace) scale(.75) legend(order(2 "Vax A" 4 "Vax B")) ti("A", position(2))

* Preference Gap *
reg gap vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined if unvaccinated == 1, nocons
coefplot, xlabel(.3 "30%" .2 "20%" .1 "10%" 0 "0%" -.1 "-10%" -.2 "-20%" -.3 "-30%" -.4 "-40%" -.5 "-50%" -.6 "-60%" -.7 "-70%") ///
ylabel(1 "No efficacy information" 2 "Symptomatic illness: 66% vs. 91%" 3 "Severe illness: 85% vs. 95%" 4 "Hospitalization/death: both > 95%" 5 "Efficacy combined") ///
graphregion(color(white)) ti("") legend(off) saving(gap_unvax.gph, replace) scale(.75) xline(0) ti("B", position(2)) xti(" " "Preference Gap (% Vax A - % Vax B)")
oneway gap group, bonferroni tab
ttest gap if efficacy_symptomatic == 1|efficacy_severe == 1, by(efficacy_severe)
ttest gap if efficacy_symptomatic == 1|efficacy_severe == 1, by(efficacy_severe)
ttest gap if efficacy_symptomatic == 1|efficacy_combined == 1, by(efficacy_combined)
oneway gap group if unvaccinated == 1, bonferroni tab
* Combined figure *
graph combine "choose_combined_unvax.gph" "gap_unvax.gph", row(2) ysize(11) xsize(8.5) graphregion(color(white))

* Figure 3 *
* Combined -- Blacks or Latino *
gen blackorlatino = 0
replace blackorlatino = 1 if black == 1
replace blackorlatino = 1 if latino == 1

reg choose_jnj vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined if blackorlatino == 1 , nocons
estimates store A
reg choose_pfizer vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined if blackorlatino == 1, nocons
estimates store B
coefplot (A, symbol(circle)) (B, symbol(square)), xlabel(0 "0%" .1 "10%" .2 "20%" .3 "30%" .4 "40%" .5 "50%" .6 "60%" .7 "70%" .8 "80%") ///
ylabel(1 "No efficacy information" 2 "Symptomatic illness: 66% vs. 91%" 3 "Severe illness: 85% vs. 95%" 4 "Hospitalization/death: both > 95%" 5 "Efficacy combined") ///
graphregion(color(white)) saving(choose_combined_nonwhites.gph, replace) scale(.75) legend(order(2 "Vax A" 4 "Vax B")) ti("A", position(2))

* Preference Gap *
reg gap vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined if blackorlatino == 1, nocons
coefplot, xlabel(.3 "30%" .2 "20%" .1 "10%" 0 "0%" -.1 "-10%" -.2 "-20%" -.3 "-30%" -.4 "-40%" -.5 "-50%" -.6 "-60%" -.7 "-70%") ///
ylabel(1 "No efficacy information" 2 "Symptomatic illness: 66% vs. 91%" 3 "Severe illness: 85% vs. 95%" 4 "Hospitalization/death: both > 95%" 5 "Efficacy combined") ///
graphregion(color(white)) ti("") legend(off) saving(gap_nonwhite.gph, replace) scale(.75) xline(0) ti("B", position(2)) xti(" " "Preference Gap (% Vax A - % Vax B)")
oneway gap group if blackorlatino == 1, bonferroni tab
ttest choose_jnj if efficacy_combined == 1, by(blackorlatino)
ttest gap if efficacy_combined == 1, by(blackorlatino)
ttest gap if efficacy_combined == 1 & (blackorlatino==1|white==1), by(blackorlatino)
ttest gap if efficacy_symptomatic == 1 & (blackorlatino==1|white==1), by(blackorlatino)
ttest gap if efficacy_severe == 1 & (blackorlatino==1|white==1), by(blackorlatino)
gen whitenotlatino = white
replace whitenotlatino = 0 if latino == 1
ttest gap if efficacy_combined == 1 & (blackorlatino==1|whitenotlatino==1), by(blackorlatino)
ttest gap if efficacy_symptomatic == 1 & (blackorlatino==1|whitenotlatino==1), by(blackorlatino)
ttest gap if efficacy_severe == 1 & (blackorlatino==1|whitenotlatino==1), by(blackorlatino)
* Combined figure *
graph combine "choose_combined_nonwhites.gph" "gap_nonwhite.gph", row(2) ysize(11) xsize(8.5) graphregion(color(white))

* Figure S1 *
* Combined -- Not Black or Latino *
reg choose_jnj vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined if blackorlatino == 0 , nocons
estimates store A
reg choose_pfizer vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined if blackorlatino == 0, nocons
estimates store B
estimates store B
coefplot (A, symbol(circle)) (B, symbol(square)), xlabel(0 "0%" .1 "10%" .2 "20%" .3 "30%" .4 "40%" .5 "50%" .6 "60%" .7 "70%" .8 "80%") ///
ylabel(1 "No efficacy information" 2 "Symptomatic illness: 66% vs. 91%" 3 "Severe illness: 85% vs. 95%" 4 "Hospitalization/death: both > 95%" 5 "Efficacy combined") ///
graphregion(color(white)) saving(choose_combined_whites.gph, replace) scale(.75) legend(order(2 "Vax A" 4 "Vax B")) ti("A", position(2))

* Preference Gap *
reg gap vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined if blackorlatino == 0, nocons
coefplot, xlabel(.3 "30%" .2 "20%" .1 "10%" 0 "0%" -.1 "-10%" -.2 "-20%" -.3 "-30%" -.4 "-40%" -.5 "-50%" -.6 "-60%" -.7 "-70%") ///
ylabel(1 "No efficacy information" 2 "Symptomatic illness: 66% vs. 91%" 3 "Severe illness: 85% vs. 95%" 4 "Hospitalization/death: both > 95%" 5 "Efficacy combined") ///
graphregion(color(white)) ti("")  saving(gap_white.gph, replace) scale(.75) xline(0) ti("B", position(2)) xti(" " "Preference Gap (% Vax A - % Vax B)")
oneway gap group if blackorlatino == 0, bonferroni tab

anova gap (group)
* Combined figure *
graph combine "choose_combined_whites.gph" "gap_white.gph", row(2) ysize(11) xsize(8.5) graphregion(color(white))

* T-tests referenced in text *
* All *
ttest choose_jnj==choose_pfizer if vax_control == 1
ttest choose_jnj==choose_pfizer if efficacy_symptomatic == 1
ttest choose_jnj==choose_pfizer if efficacy_severe== 1
ttest choose_jnj==choose_pfizer if efficacy_hospitalization == 1
ttest choose_jnj==choose_pfizer if efficacy_combined == 1

* Unvaccinated *
ttest choose_jnj==choose_pfizer if vax_control == 1 & unvaccinated == 1
ttest choose_jnj==choose_pfizer if efficacy_symptomatic == 1 & unvaccinated == 1
ttest choose_jnj==choose_pfizer if efficacy_severe== 1 & unvaccinated == 1
ttest choose_jnj==choose_pfizer if efficacy_hospitalization == 1 & unvaccinated == 1
ttest choose_jnj==choose_pfizer if efficacy_combined == 1 & unvaccinated == 1

* Whites *
ttest choose_jnj==choose_pfizer if vax_control == 1 & white == 1
ttest choose_jnj==choose_pfizer if efficacy_symptomatic == 1 & white == 1
ttest choose_jnj==choose_pfizer if efficacy_severe== 1 & white == 1
ttest choose_jnj==choose_pfizer if efficacy_hospitalization == 1 & white == 1
ttest choose_jnj==choose_pfizer if efficacy_combined == 1 & white == 1

* Non-Whites *
ttest choose_jnj==choose_pfizer if vax_control == 1 & blackorlatino== 1
ttest choose_jnj==choose_pfizer if efficacy_symptomatic == 1 & blackorlatino== 1
ttest choose_jnj==choose_pfizer if efficacy_severe== 1 & blackorlatino== 1
ttest choose_jnj==choose_pfizer if efficacy_hospitalization == 1 & blackorlatino== 1
ttest choose_jnj==choose_pfizer if efficacy_combined == 1 & blackorlatino== 1

ttest choose_jnj if efficacy_symptomatic == 1 & blackorlatino== 1 | efficacy_combined == 1 & blackorlatino== 1, by(efficacy_combined)
ttest choose_jnj if efficacy_severe == 1 & blackorlatino== 1 | efficacy_combined == 1 & blackorlatino== 1, by(efficacy_combined)
ttest choose_jnj if efficacy_combined == 1, by(white)

* Preference Gaps *
gen gap = choose_jnj-choose_pfizer
* All respondents; SI Figure 1 *
reg gap vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined , nocons
coefplot, xlabel(.2 "20%" .1 "10%" 0 "0%" -.1 "-10%" -.2 "-20%" -.3 "-30%" -.4 "-40%" -.5 "-50%" -.6 "-60%" -.7 "-70%") ///
ylabel(1 "No efficacy information" 2 "Symptomatic illness: 66% vs. 91%" 3 "Severe illness: 85% vs. 95%" 4 "Hospitalization/death: both > 95%" 5 "Efficacy combined") ///
graphregion(color(white)) ti("") legend(off) saving(gap_all.gph, replace) scale(.75) xline(0)
ttest gap if efficacy_symptomatic == 1|efficacy_severe == 1, by(efficacy_severe)
ttest gap if efficacy_symptomatic == 1|efficacy_combined == 1, by(efficacy_combined)

* Non-whites *

reg choose_jnj vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined if white == 0 , nocons
estimates store A
reg choose_pfizer vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined if white == 0, nocons
estimates store B
coefplot (A, symbol(circle)) (B, symbol(square)), xlabel(0 "0%" .1 "10%" .2 "20%" .3 "30%" .4 "40%" .5 "50%" .6 "60%" .7 "70%" .8 "80%") ///
ylabel(1 "No efficacy information" 2 "Symptomatic illness: 66% vs. 91%" 3 "Severe illness: 85% vs. 95%" 4 "Hospitalization/death: both > 95%" 5 "Efficacy combined") ///
graphregion(color(white)) saving(choose_combined_nonwhites.gph, replace) scale(.75) legend(order(2 "Vax A" 4 "Vax B")) ti("A", position(2))

* Preference Gap *
reg gap vax_control efficacy_symptomatic efficacy_severe efficacy_hospitalization efficacy_combined if white == 0, nocons
coefplot, xlabel(.3 "30%" .2 "20%" .1 "10%" 0 "0%" -.1 "-10%" -.2 "-20%" -.3 "-30%" -.4 "-40%" -.5 "-50%" -.6 "-60%" -.7 "-70%") ///
ylabel(1 "No efficacy information" 2 "Symptomatic illness: 66% vs. 91%" 3 "Severe illness: 85% vs. 95%" 4 "Hospitalization/death: both > 95%" 5 "Efficacy combined") ///
graphregion(color(white)) ti("") legend(off) saving(gap_nonwhite.gph, replace) scale(.75) xline(0) ti("B", position(2)) xti(" " "Preference Gap (% Vax A - % Vax B)")
oneway gap group if blackorlatino == 1, bonferroni tab
ttest choose_jnj if efficacy_combined == 1, by(blackorlatino)
oneway gap group if white == 0, bonferroni tab

* Combined figure *
graph combine "choose_combined_nonwhites.gph" "gap_nonwhite.gph", row(2) ysize(11) xsize(8.5) graphregion(color(white))

